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Abstract 

We simulate single and multiple Ising models coupled to 2~d gravity 
using both the Swendsen-Wang and Wolff algorithms to update the 
spins. We study the integrated autocorrelation time and find that 
there is considerable critical slowing down, particularly in the mag- 
netization. We argue that this is primarily due to the local nature 
of the dynamical triangulation algorithm and to the generation of a 
distribution of baby universes which inhibits cluster growth. 
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Considerable work has been devoted to the study of the pgrlormance 
of cluster algorithms in reducing the critical slowing down of niaiy statis- 
tical models. The Swendsen-Wang [1] and the Wolff [2] algoritfms have 
proven very effective in beating the critical slowing down (CSDjexhibited 
by these models when simulated with a standard Metropolis [3] update. 
Reviews on cluster algorithms and CSD may be found on refs. [4]. In 
recent years, the numerical study of Potts models coupled to 2-d gravity 
has received much attention [5, 6, 7], aided as well by increased analytical 
understanding of these models [8] . Cluster algorithms have proved useful 
in saving computational effort on the update of the Potts variables. The 
present study is motivated by the fact that there is little understanding 
of the actual extent of the improvement achieved in these simulations. It 
is worthwhile, therefore, to measure CSD in the case of simple Ising spins 
coupled to a dynamical lattice in order to quantify the performance of 
cluster algorithms. We find that there is considerable CSD, especially in 
the magnetization, and we relate this to the dynamics of cluster formation 
on a random lattice. 



£661 AON £Z ' IA9eo"l ie6/l^I-"d3M-"AIXJ^: 



We shall consider a model in which Ug Ising spins are attached to the 
vertices of triangulations. The triangulations are characterized by their 
adjacency matrix Cy, which equals 1 if i and j are neighbors and vanishes 
otherwise. Cy is the discrete analogue of the world-sheet metric Qiy We 
shall restrict ourselves to the set of triangulations with vertices T/v 
containing only loops of length 3 or greater and vertices of coordination 
number of at least 3. The triangulation has a fixed toroidal topology. We 
simulate a theory determined by the partition function 
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where a labels the spin species. In refs. [9, 10] we investigated in detail 
this model for the cases Ug = 1 and 2. We measured spin susceptibility and 
percolation critical exponents using finite-size scaling and showed that 
logarithmic corrections to scaling were essential for agreement between 
the measured and theoretical exponents. In this paper, we deal with the 
cases Tig = 1, 2 and 3, concentrating on the issue of CSD and its origin. 

The standard way of implementing the partition function (1) via a 
Monte Carlo simulation is to use the Swendsen-Wang (SW) or Wolff clus- 
ter algorithm to update the spin variables and to use the "link flip" [11] 
to simulate the sum over all triangulations. To implement a SW update 
one first divides all of the spins into Fortuin-Kasteleyn (FK) clusters [12]. 
These clusters of bonded spins are created by introducing bonds between 
same sign spins with probability p = 1 — exp(— 2/3). Then, one flips all 
clusters with probability one-half. The Wolff algorithm consists of ran- 
domly choosing a spin, constructing a FK cluster around it and flipping 
the cluster with probability one. To compare the autocorrelation times 
of the SW and Wolff algorithms one should define the Wolff update so 
that both algorithms require comparable CPU time. For this reason we 
choose a Wolff update to consist of consecutive flips of FK clusters that 
reverse the sign of at least 40% of the spins. One alternative to this would 
be to scale the correlation times using the average cluster size [13]. Each 
spin update precedes a mesh update, in which we attempt to flip 3N ran- 
domly chosen links, being the number of vertices of the triangulation. 
Our implementation ensures that the relative number of mesh and spin 



updates is roughly equivalent. 

The observables that we analyzed are the energy density, the mag- 
netization density, the susceptibility (namely the magnetization squared) 
and the average value of |g — 6|, where q is the coordination number of 
a vertex of the triangulation. In some of the simulations that employed 
SW updates, we also measured the mean size of pure percolation and FK 
clusters. The mean FK cluster size Svk is given by the quantities (s)woifr 
and (s^)sw/('5)sw; s denotes the number of sites constituting a cluster 
and averages are taken over the distribution of clusters built in the Wolff 
and SW algorithms respectively. The magnetic Ising observables are di- 
rectly related to the structure of these clusters [14]. In particular, for 
/5 < /3c, 'S'fk is equal to the susceptibility, defined as 

X = |(M2), □ □ (2) 

where M denote5^?he average magnetization density. Ac^ujilly, this alter- 
native definition of x is used as a reduced variance estimator [15]. 

Before moving on to present our results, we discuss how we esti- 
mated the autocorrelation times. It is known [16] that the following 
relation holds between the estimators of the naive and true variance of 
an observable O; 

^ Var(0)truc ^ 2 Tint Var(0)naive- (3) 

For a subset of our data, we also directly measured the autocorrelation 
function and through standard methods [16] computed Tint. We verified 
(see also [17]) that both of the above techniques gave consistent values of 
Tint. The variances of the observables were extracted from the raw data 
using the binning method. To extract the dynamic exponent, we fit the 
data according to the scaling ansatz 
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It is difficult to estimate the linear size of random triangulations. There- 
fore, we shall leave explicit the dependence of the exponent upon the 
intrinsic Hausdorff dimension dn of the triangulation and fit our data 



2 



using the total area N. The values of z/dfj were extracted from the 
auto-correlation data using a log- log regression fit, excluding from the fit 
the results for the smaller volumes, since they are affected the most by 
finite-size effects. The values of P^r degree of freedom were always of 
order one. In the case of the magnetization for the = 2 and 3 models, 
Wolff algorithm, we extracted the exponent using only the two largest 
lattices, since there are larger finite-size effects. 





z/du 






Model 


Metropolis 


SW 


Wolff 


ns = 1 


.85 ± .06 


.58 ± .05 


.54 ± .05 


ns = 2(*) 


.95 ± .05 


.62 ± .06 


.58 ± .09 


ns = 3(*) 


.9± .1 


.49 ± .08 


.55 ± .1 



Table 1: Critical exponent z/dn for the Magnetization from fits. 
(*) These numbers are not reliable; we discuss this point in the text. 

We present results for three different models, ng = 1, 2 and 3. Each 
model was simulated with the SW, Wolff and Metropolis algorithms. The 
ng = 1 model was simulated at the critical value of (3, which is known an- 
alytically [9] . For the other models we chose the /3 value by looking at the 
peak of the susceptibility and the intersection of the Binder's cumulant 
[18] curves [9]. Each simulation (model and algorithm) was run at four or 
five values of N (512, 1024, 2048, 4096 and 8192 in the ns = 1 model) and 
each consisted of 10^ thermalization sweeps and 3-5x10^ measurement 
sweeps. Measurements were taken every sweep. 

Prom the analysis presented in figures 1, 2, 3 and 4 and tables 1 and 
2, we deduce the following: 

1. There is considerable critical slowing down in these models. Figure 1 
shows the improvement gained by the use of cluster algorithms — the 
autocorrelation times and the dynamic exponent are significantly 
lower than the corresponding Metropolis values. 

2. In all cases, the magnetization is the observable that suffers most 
from critical slowing down. This behavior is quite different from 
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Model 


Metropolis 


SW 


Wolff 


"-8=1 


.62 ± .03 


.057 ± .005 


.04 ± .03 


"-8 = 2 


.35 ± .1 


.08 ± .02 


.17± .08 


"s = 3 


.5±.l 


.05 ± .04 


.37± .08 



Table 2: Critical exponent z/dn for the Energy from fits. 

that of the Ising model on a regular lattice, where the energy exhibits 
CSD equal to or greater than that of the magnetization [19]. In our 
simulations the observable |g — 6| did not show any significant CSD. 

The SW and Wolff algorithms have similar 'performance within the 
statistical accuracy of our data. Qfi-|the smaller lattices the Wolff 
algorithm is somewhat more efficient than SW, but this advantage is 
a finite size artifact since it disappears on the larger lattices. For 2-d 
Ising models on flat and Poissonian lattices these algorithms exhibit 
roughly comparable performance [17, 19]. In the 3~d case, the Wolff 
algorithm is much more efficient. iReems that the relative perfor- 
mance of SW and Wolff algorithms depends on the dimensionality of 
the lattice. !!□ 

It is hard to determine differences iizthe degree of CSD between one 
and two Ising models coupled to gravity, given our statistics. In [9] 
we found, in fact, that the numerically measured behavior of the one 
and two species models is very similar. This, however, is a conse- 
quence of logarithmic corrections [20] in the two-species model. One 
might suspect that the scaling law for r in this model incorporates 
logarithms as well. On much larger lattices where logarithmic be- 
havior is distinguishable from small power law scaling, the effective 
CSD might be considerably different in the two generation, as com- 
pared to one generation, case. This situation is similar to that of the 
2-d 4-state Potts model. Here, Li and Sokal [21] have shown that 
measurements of ^ obtsilied using the ansatz (4) , violate rigorous 
bounds. They suggest that these measurements of z are not correct 
because the fits to r fail to take into account logarithmic corrections. 
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Figure 1: Comparison of the integrated autocorrelation times for the mag- 
netization in the Ug = 1 model. The dashed lines are log-log regression 
fits. 

For the ng = 3 model the situation is worse; the corrections to scal- 
ing of (4) may be even larger. In this case, there are no theoretical 
arguments that predict the form of these corrections. Therefore, we 
anticipate that the numbers we have quoted for z/dn differ consid- 
erably from the correct asymptotic values in the ng = 2 and 3 cases. 
To give a sense of the magnitude of the corrections to scaling we note 
that on similar size lattices, the estimate of ^/vdn differs by about 
3% from its asymptotic value in the ng = 1 case and by almost 50% 
in the ng = 2 case [9]. 

The presence of CSD in these models should not be a surprise, since 
the triangulation is updated locally. In this context, it is relevant to 
briefly recall the results of a similar analysis [10] in which percolation 
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Figure 2: Comparison of the integrated autocorrelation times for the 
energy in the rig = 1 model. 

□ 

clusters were studied on random triangulations without matter (pure 
gravity). The meshes were updated locally; all link flips that did not 
lead to degenerate triangulations were allowed. It has been shown that 
these random triangulations are characterized by a scaling distribution 
of baby universes [22]. The formation of clusters is quite sensitive to 
the presence of the bottlenecks (see fig. 5), which inhibit cluster growth 
into and out of baby universes. Since this structure of baby universes is 
slow to decorrelate under the local link-flip updatef^^e expect that the 
mean percolation cluster sizes will be afflicted by critical slowing down. 
In the case of pure gravity, the mean size of percolation clusters built 
on these triangulations in fact exhibits critical slowing down of magni- 
tude zjdn = .70(2). We also cojEsSructed percolation clusters in the Ising 
simulations that used SW spin updates. We observed, in this case, a 
similar CSD of the world-sheet geometry; zjdu for the mean size of pure 
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Figure 3: Comparison of the integrated autocorrelation times for the 
magnetization in the ng = 2 model. We excluded smaller volume points 
from the fits. 



percolation clusters was measured to be .74(6). 

We now argue that the critical slowing down in the gravity sector 
should lead to considerable CSD for magnetic observables. For one would 
expect that baby universes should trap FK clusters as well as pure per- 
colation clusters. As said before, for an Ising model simulated on an 
arbitrary random triangulation, the mean FK cluster size equals the av- 
erage magnetic susceptibility. The value of the magnetization is thus 
clearly sensitive to those features of the geometry that strongly affect 
the FK cluster size. This coupling transfers critical slowing down to the 
magnetic sector. Some evidence in support of this argument follows from 
our measurements of the mean size of the FK clusters built to perform 
SW updates. This observable exhibited a value of z/dn of .52(6). 
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Figure 4: Comparison of the integrated autocorrelation times for the 
energy in the ng = 2 model. 

The efficiency of cluster algorithms is also typically affected by the 
distribution of cluster sizes. If the clusters are too small, ffipping them 
will fail to decorrelate distant spins. If one cluster fills most of the lattice, 
successive flips will essentially cancel each other out. Indeed FK clusters, 
on average, are miich smaller in the dynamical case than in the case of a 
fixed flat lattice. Their mean size is determined by the exponent ^/vdn', 
>Sfk ~ ATT/'^rfH. For flat lattices, -^/udu = 7/8. It is 2/3 for the ng = 1 
Ising model, which is quite close to the value for the 3-d Ising model^. 

It is difficult to determine directly whether the autocorrelation times 
are influenced primarily by the slow decorrelation of the world-sheet ge- 
ometry or the effects of the smaller cluster sizes. The argument, however, 

^Simulations of the 3-d Ising model with cluster algorithms suffer noticeable critical 
slowing down, though still rather modest compared to the CSD observed here. In the 
flat 3-d case the energy is the observable showing the greater CSD. 




Figure 5: Schematic representation of a baby universe. 



that CSD arises from updates of the geometry apphes only to magnetic 
observables, and not to the energy density. Since the CSD we observe is 
greater for the magnetization, it seems hkely that the slow decorrelation 
of the world-sheet geometry is the primary mechanism responsible for 
CSD. 
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